Compare urban traits and phylogency

source('../env.R')

Using GitHub PAT from the git credential store.
Skipping install of 'clootl' from a github remote, the SHA1 (30ac57ee) has not changed since last install.
  Use `force = TRUE` to force installation

Create dataset

Load community data

communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))

Rows: 2428 Columns: 7── Column specification ────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, ebird_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

communities

Load trait data

traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_ebird.csv'))

Rows: 332 Columns: 10── Column specification ────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): ebird_species_name, habitat, trophic_level, trophic_niche, primary_lifestyle
dbl (5): beak_width, hwi, mass, habitat_density, migration
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(traits)

species_urban_response = communities %>% group_by(ebird_species_name) %>% summarise(mean_relative_abundance_proxy = mean(relative_abundance_proxy)) %>% left_join(traits)

Joining with `by = join_by(ebird_species_name)`

species_urban_response$is_urban = species_urban_response$mean_relative_abundance_proxy > 0
species_urban_response$ebird_species_name = str_replace(species_urban_response$ebird_species_name,' ', '_')
species_urban_response

Most urban species

species_urban_response %>% filter(mean_relative_abundance_proxy > 50)

species_urban_response %>% filter(mean_relative_abundance_proxy > 40)

species_urban_response %>% filter(mean_relative_abundance_proxy > 30)

Difference amongst urban and non-urban

species_urban_response %>% group_by(is_urban) %>% summarise(
  median_beak_width = median(beak_width),
  iqr_beak_width = IQR(beak_width),
  
  median_hwi = median(hwi),
  iqr_hwi = IQR(hwi),
  
  median_mass = median(mass),
  iqr_mass = IQR(mass),
  
  n = n()
)

Compare traits between urban and non-urban species

Gape Width

ggplot(species_urban_response, aes(x = is_urban, y = beak_width)) + geom_boxplot()

wilcox.test(beak_width ~ is_urban, species_urban_response)

    Wilcoxon rank sum test with continuity correction

data:  beak_width by is_urban
W = 1525, p-value = 0.3893
alternative hypothesis: true location shift is not equal to 0

There is not a significant response for the gape width between urban species and non urban species.

HWI

ggplot(species_urban_response, aes(x = is_urban, y = hwi)) + geom_boxplot()

wilcox.test(hwi ~ is_urban, species_urban_response)

    Wilcoxon rank sum test with continuity correction

data:  hwi by is_urban
W = 1213, p-value = 0.01492
alternative hypothesis: true location shift is not equal to 0

There is a significant response between urban and non-urban species

Mass

ggplot(species_urban_response, aes(x = is_urban, y = mass)) + geom_boxplot()

wilcox.test(mass ~ is_urban, species_urban_response)

    Wilcoxon rank sum test with continuity correction

data:  mass by is_urban
W = 1825, p-value = 0.5168
alternative hypothesis: true location shift is not equal to 0

There is not a significant response between species.

Phylogeny

Phylogentic response to relative abundance?

tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
tree_pruned = ladderize(drop.tip(tree, setdiff(tree$tip.label, species_urban_response$ebird_species_name)))
ggtree(tree_pruned, layout="fan")

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

columbidae_response <- species_urban_response$mean_relative_abundance_proxy
names(columbidae_response) <- species_urban_response$ebird_species_name

Is the a phylogenetic signal for our relative abundance proxy

col_ut_phylo_signal <- phylosig(tree_pruned, columbidae_response, method="lambda", test=TRUE)
col_ut_phylo_signal

Phylogenetic signal lambda : 0.46832 
logL(lambda) : -576.839 
LR(lambda=0) : 7.45324 
P-value (based on LR test) : 0.00633224 

Plot tree

species_urban_response$ebird_species_name_clean = gsub("_", " ", species_urban_response$ebird_species_name)

normalize <- function(x) (x - min(x)) / (max(x) - min(x)) * 5 + 1  # Scale to range [1, 6]
species_urban_response$beak_width_norm <- normalize(species_urban_response$beak_width)
species_urban_response$hwi_norm <- normalize(species_urban_response$hwi)
species_urban_response$mass_norm <- normalize(species_urban_response$mass)

g = ggtree(tree_pruned, layout="fan") %<+% species_urban_response + abundance_proxy_scale

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

g + 
  geom_fruit(geom = geom_point, mapping = aes(size = beak_width_norm), fill = beak_width_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0, show.legend = FALSE) + 
  geom_fruit(geom = geom_point, mapping = aes(size = hwi_norm), fill = hwi_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0.13, show.legend = FALSE) + 
  geom_fruit(geom = geom_point, mapping = aes(size = mass_norm), fill = mass_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0.15, show.legend = FALSE) + 
  geom_tiplab(size=2, aes(label=ebird_species_name_clean, color = mean_relative_abundance_proxy), offset = 15) +
  labs(color = "Mean relative abundance proxy") + 
  theme(legend.position = "bottom") + ggplot2::xlim(0, 80)

ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny.jpg'), width = 2000, height = 2100, units = 'px')

Format plot for nature

g_fig_base = ggtree(tree_pruned, layout="fan") %<+% species_urban_response

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

g_fig = g_fig_base + 
  geom_fruit(geom = geom_point, mapping = aes(size = beak_width_norm), fill = beak_width_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0, show.legend = FALSE) + 
  geom_fruit(geom = geom_point, mapping = aes(size = hwi_norm), fill = hwi_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0.13, show.legend = FALSE) + 
  geom_fruit(geom = geom_point, mapping = aes(size = mass_norm), fill = mass_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0.15, show.legend = FALSE) + 
  geom_tiplab(size=2, aes(label=ebird_species_name_clean), offset = 30)

heatmap_fig_ata = data.frame(species_urban_response[, "mean_relative_abundance_proxy"])
colnames(heatmap_fig_ata) = c('MRAP')
rownames(heatmap_fig_ata) = species_urban_response$ebird_species_name

gheatmap(g_fig, 
         heatmap_fig_ata, 
         offset = 14, 
         width = 0.2,
         colnames_angle = 0, 
         colnames_offset_x = 0.5, 
         font.size = 3) + 
  abundance_proxy_scale_fill + 
  labs(fill = "MRAP (Mean relative abundance proxy)") + 
  theme(legend.position = "bottom") + 
  ggplot2::xlim(0, 85)

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

ggsave(filename(FIGURES_OUTPUT_DIR, 'figure3.jpg'), width = 180, height = 200, units = 'mm', dpi = 450)

Plot traits against abundance proxy

cut_off_3 = 50
cut_off_2 = 35
cut_off_1 = 1

hwi_by_mass_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(hwi, mass))
hwi_by_mass_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(hwi, mass))
hwi_by_mass_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(hwi, mass))

beak_width_by_mass_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(beak_width, mass))
beak_width_by_mass_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(beak_width, mass))
beak_width_by_mass_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(beak_width, mass))

hwi_by_beak_width_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(hwi, beak_width))
hwi_by_beak_width_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(hwi, beak_width))
hwi_by_beak_width_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(hwi, beak_width))

species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% nrow()
species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% nrow()
species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% nrow()

polygon_line_type = 'dashed'
polygon_linewidth = 0.4
polygon_alpha1 = 0.05
polygon_alpha2 = 0.1
polygon_alpha3 = 0.15

hwi_by_mass = ggplot(species_urban_response, aes(x = hwi, y = mass, colour = mean_relative_abundance_proxy)) + 
  geom_polygon(data = hwi_by_mass_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = hwi_by_mass_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = hwi_by_mass_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_point() + 
  abundance_proxy_scale +
  labs(colour = "Mean relative abundance proxy", y = 'Mass', x = '') + 
  theme(legend.position = 'bottom', legend.title.position = 'top', legend.key.width = unit(12, 'mm'))

beak_width_by_mass = ggplot(species_urban_response, aes(x = beak_width, y = mass, colour = mean_relative_abundance_proxy)) + 
  geom_polygon(data = beak_width_by_mass_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = beak_width_by_mass_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = beak_width_by_mass_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_point() + 
  abundance_proxy_scale + 
  labs(y = '', x = 'Beak width')

hwi_by_beak_width = ggplot(species_urban_response, aes(x = hwi, y = beak_width, colour = mean_relative_abundance_proxy)) + 
  geom_polygon(data = hwi_by_beak_width_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = hwi_by_beak_width_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = hwi_by_beak_width_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_point() + 
  abundance_proxy_scale + 
  labs(y = 'Beak width', x = 'HWI')

legend <- ggpubr::get_legend(
  # create some space to the left of the legend
  hwi_by_mass
)

plot_grid(
  nrow = 2, ncol = 2,
  hwi_by_mass + theme(legend.position="none"),
  beak_width_by_mass + theme(legend.position="none"),
  hwi_by_beak_width + theme(legend.position="none"),
  legend
)

ggsave(filename(FIGURES_OUTPUT_DIR, 'traits.jpg'), width = 2000, height = 1800, units = 'px')